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Abstract 

We consider hyperbolic systems of conservation laws which are discretized in space 
by spectral collocation methods and advanced in time by finite difference schemes. At 
any time-level we introduce a domain decomposition method based on an iteration-by- 
subdomain procedure yielding at each step a sequence of independent subproblems (one 
for each subdomain) that can be solved simultaneously. The method is set for a general 
nonlinear problem in several space variables. The convergence analysis, however, is 
carried out only for a linear one-dimensional system with continuous solutions. A 
precise form of the error reduction factor at each iteration is derived. Although the 
method is applied here to the case of spectral collocation approximation only, the 
idea is fairly general and can be used in a different context as well. For instance, its 
application to space discretization by finite differences is straightforward. 


1 Research was suppported under the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS 1-18605 while the author was in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 
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Introduction 

In this paper we propose an effective domain decomposition method for the numerical solution of 
quasi-linear hyperbolic systems of conservation laws. 

Though our emphasys is on spatial discretization by spectral collocation schemes, the domain de- 
composition approach we advocate here can be adopted for finite difference discretizations as well. 
At every time-level the method applies after that the initial value problem has been advanced by 
either an implicit or an explicit time-stepping. 

The spatial computational domain is subdivided into several adjoining, non intersecting subdo- 
mains; within each of them we look for a polynomial solution of degree N with respect to each 
component. If T, denotes the interface between two of such subdomains, say Q, and Q,+i , at each 
gridpoint on T,- we enforce the conditions of continuity of the physical variables. Furthermore we 
require the fulfillment of the so-called compatibility equations, i.e., of those characteristic combi- 
nations of the original equations which express wave propagation across T, from Q, toward Q,+i 
and viceversa. 

The main reason for using the multidomain spectral method rather than the standard single-domain 
spectral approach stems from its capability of covering problems in complex geometry. More- 
over, the spectral multidomain approach allows local refinement to resolve internal layers (or even 
discontinuities) maintaining however the spectral accuracy enjoyed by the classical spectral collo- 
cation method (e.g. [CHQZ], Ch.12). 

We also propose here an iteration-by-subdomain algorithm which allows the decoupling of the 
subproblems arising from the multidomain approach making it possible to solve at each iteration 
as many independent subproblems as the number of subdomains. This algorithm requires that on 
T, the values of the characteristic variables impinging Q, equate those outgoing from Q, + i at the 
previous iteration, and viceversa. 

For the case of linear hyperbolic systems we prove that the above iteration-by- subdomain method 
is convergent. Furthermore, we find a close expression of the reduction factor per iteration and we 
show that it is independent of N, the number of gridpoints inside each subdomain. Furthermore, an 
algebraic interpretation of our iteration-by-subdomain algorithm is derived in terms of the Schur 
complement (or capacitance) matrix. 

Finally, we show how the method proposed here can be adapted to the case in which an internal 
interface is a discontinuity surface (actually, a "k-shock" in the sense of Lax [L]). This typically 
occurs if a shock fitting approach is adopted ([CHQZ], sect.8.6). In such a case, the compatibility 
equations together with the Rankine-Hugoniot conditions allow us to compute first the flow field 
within the upstream subdomain then the one in the downstream subdomain. 

Examples of domain decomposition methods of similar type have been proposed in the latest years. 
We recall for instance [K] and the references given in [CHQZ], Ch.13. We also refer to [MS] and 
[MSH] for applications of spectral multi-domain techniques to reacting flow whose shape and mo- 
tion are generated at each time-level. In the frame of finite differences we refer to the earlier paper 
[CGV V] where the issue of the compatibility equations at subdomain interfaces was adressed. More 
recent applications of finite difference multidomain methods for compressible flow simulations can 
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be found in the reference [CGPW]. 

An outline of this paper is reported below. In section 1 we introduce the initial-boundary-value 
problem and the compatibility equations, while in section 2 we state the domain decomposition 
formulation of the Chebyshev collocation approximation to the problem. In section 3 we present an 
iteration-by-subdomain algorithm for an effective solution of the domain decomposition problem. 
Thereafter we confine to the case of one-dimensional hyperbolic systems. We write the domain 
decomposition problem as well as the iteration-by-subdomain algorithms. In section 5 we carry 
out the convergence analysis for the iterative algorithm, and in section 6 we derive the relationship 
with the Schur complement matrix associated with the interface unknowns. Finally, in section 7, 
we report some numerical results for a one-dimensional model problem. 

1. The Hyperbolic system and the compatibility equations 

2. Spectral collocation approximation to the system (1.1): a domain decomposition approach 

3. An iteration by subdomain algorithm for the solution of the domain decomposition problem 

4. A particular case: one dimensional problems 

5. Convergence analysis for the spectral collocation problem 

6. Capacitance matrix interpretation 

7. Numerical Results. 


1. The Hyperbolic System and the Compatibility Equations 
We consider the following differential system 


( 1 . 1 ) 


du ^ A , du , . 
gl + 22 A i( u ) — = f » n 
y=i 


dxj 


Q x (0, T) 


where m = 1, 2 or 3, Q is an open bounded domain of R m , u and / are two vector functions: 
u, / : £2 x (0, T) — ♦ R p , with p > 1, and Aj(u) are p x p matrices possibly depending on u. 

The system (1.1) is supplemented by initial conditions of the form: 


(1.2) u(x, 0) = p(x) Vx&Cl 

and by suitable boundary conditions at (0, T) x dQ that will be specified later on. 

For any £ E R m such that |£| = £/) = 1* let us define the characteristic matrix for the £ 

direction as follows: 


M f) := Y, A i (u> b 

;=i 


(1.3) 
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The system (1.1) is assumed to be hyperbolic in time: this means that for any such A(() has p 
real eigenvalues and moreover it is diagonalisible. 

Let us denote by {A*, k = 1, ...,p} the eigenvalues of A(0» and by (v fe , k = 1, ..., p} the set of the 
corresponding left eigenvectors, so that 


(1.4) v k A(£) = A k v k fc=l,...,p 

We assume that A* > 0 for k - 1,...,? and A* < 0 for k = q + 1, ...,p for some 0 < q < p 
(obviously, q depends on the particular direction £ that we are considering). Let us take the inner 
product of v k with (1.1); for each k we obtain the scalar equation 


(1.5) 



k- 1.....P 


Denote now by t the direction orthogonal to so that r • £ = 0, |r| = 1. Since 


(L6) 


du du du 
dxj d£ * + dr T * 


j = 1, ...,m 


from (1.5) and (1.4) it follows that: 


(1.7) 


v" 


(du 



Aj(u) 


du 

dr 



k = 1 , 


P 


Following a terminology proposed in [CGVV], we will refer to (4.7) as to the compatibility equa- 
tions for the problem (1.1). 


Remark 1.1 Let T denote a (m-l)-dimensional manifold of R m , and denote by v the unit normal 
direction to T. Taking £ = v, equations (1.7) restricted to T become: 


( 1 . 8 ) 




fc= l,...,p 


The right hand side of (1.8) depends on the tangential derivatives of u on T, while the left hand side 
yields a combination (through the components of the eigenvector v k ) of transport equations along 
a direction which is normal to T. 
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Assume T is the boundary of £2 and v is oriented outward 12 (see fig. 1.1 a); then for A: = 1, q, 
(1.8) yield q transport equations according to which information are propagated from the inside to 
the outside of 12. For such a reason, equations (1.8) for k = 1, ..., q will be called the compatibility 
equations for the domain 12 . 

If we now assume that 12i and 122 are two adjoining subdomains of 12, and T is their common 
boundary, taking as u the normal direction to T, oriented from 12i , to Q 2 (see fig. 1 . 1 b), then the first 
q equations of (1.8) are the compatibility equations for 12i (they entail propagation of information 
from 12i to 122). Obviously, for k = q + 1, ...,p, (1.8) provide the compatibility equations for 12? . 
Such a distinction will be crucial in the definition of the multidomain method we are going to 
propose in the forthcoming sections.o 



Fig. 1.1 a 



2. Spectral Collocation Approximation to the System (1.1): a Domain Decomposition Ap- 
proach 

We assume in this section that 12 is an m-dimensional hypercube. Let 12i and CI 2 be two open 
subregions so that 12 = 12 i U 122, 12 i fl 122 - T, 12 i D 122 = < f > where T is orthogonal to one cartesian 
direction. We assume here that the solution of (1.1) is continuous across T. The case of solutions 
that are discontinuous across T is faced in sect.2.1. 

Let N be a given positive integer; we denote by 


(2.1) E := |atk = cos— , k = (ki,...,k m ) , 0 < fc, < N , » = 

the set of Chebyshev-Lobatto points of the reference hypercube [-1, l] m . Further we denote by E 1 
and £ 2 the corresponding set of points in 12i and Q 2 respectively. Note that a point of E 1 (resp.E 2 ) 
lies on the boundary of 12i (resp.122) if at least one of its subindices kj is either 0 or N. Since we 
are using the same number of points in 12i and Q 2 . the points of E 1 U T and those of E 2 U T are 
exacdy the same. We will denote by Ep this common set of interface points. 
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A (continuous in time) spectral domain decomposition method for problem (1.1) is defined as 
follows. At each time t > 0 we look for u' N (t) e (PN(fli)) p that satisfies the following set of 
equations. 

(a) At each point of X 1 internal to Qj 


( 2 . 2 ) 




y=i 


(b) at each point of X 2 internal to Clz 


(2.3) 

(c) at each point of Xp 




dt 


;=i 


(2.4) «ir = »5r 

Moreover we require that 



v is the outward unit normal direction to Qi on T (see fig. 1.1 b). A* and v k are the eigenvalues and 
the (left) eigenvectors of the matrix A{y) (with A* > 0 for k = 1, ..., q), and r is the tangential 
direction on T. 

Furthermore we satisfy the equations: 


( 2 . 6 ) 


v 




du% 
i=i 


dr 


■) 


k = q+l,...,p 


Following the definitions given in the Remark 1.1, (2.4) are the q compatibility equations for fii 
on r, while, by virtue of (2.4), (2.5) are the p — q ones for 


(d) at each point of X 1 belonging to a "face" 0 of (excluding the interface T) we set 
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(2.7) 



y=i 



fc = i, 


9 


where i/ is now the outward normal direction to Qi on O, r is the tangential direction on d>, while 
A* and v k are the eigenvalues and the (left) eigenvectors of the matrix A(u) (note that here q is not 
necessarily the same of (2.5)). 

The remaining p — q equations that are needed at each collocation point of O must be provided by 
the physical boundary conditions that supplement (1.1) and (1.2). 


(e) at each point of E 2 belonging to a "face" O of <9Q 2 (excluding the interface T) we enforce 


( 2 . 8 ) 


v 



du N 

at 


+ A 


du 


*) 


V” 


| /- $>,(«&) 


du N 

dr 



k = i, 
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where v is the outward normal direction to 12 2 on C>, t the tangential direction on T, X k and v k are 
the eigenvalues and the (left) eigenvectors of the matrix A(v). (Again, X k > 0 for k = 1, ...,q, 
where q here is not necessarily the same as before). 

The remaining p — q equations are prescribed as boundary conditions. 


Remark 2.1 In the frame of classical single-domain spectral collocation methods for hyperbolic 
systems, the use of compatibility equations for collocation boundary points was firstly advocated in 
[GGT] and in [CQ]. For the same method, a stability and convergence analysis were developed in 
[GLT1 ] and [GLT2] for the case of dissipative boundary conditions. No convergence analysis (with 
respect to N) is available so far for spectral multidomain approximations of hyperbolic systems.^ 


2.1 The case of discontinuous solutions across the subdomain interface 
We consider now the case in which T not an arbitrary surface anymore, but rather it is precisely the 
front shock propagating throughout the computational domain Q. This case typically occurs when 
a shock fitting technique is adopted with the purpose of determining the shape and the motion law 
of the shock front at each time-level (see, e.g. [CHQZ], sect.8.6). We still denote by Qi and Q 2 the 
adjoining subdomains separated by T, by v the outward normal vector to Qi on T and with w the 
speed with which T propagates in the direction u. 

We will assume that there exists a vector function F: Qx (0, T) — » 1R P such that Aj(u) = dF ^ - 
so that (1.1) can be written in conservation form as follows: 


du A 8Fj(m) 

Q t JL, du 1 


in Q x (0, T) 


(2.9) 
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The weak solutions to (2.9) (e.g., [L]) are allowed to be discontinuous across the interface T, but 
should satisfy the following jump conditions ( Rankine-Hugoniot conditions ) 


m 

(2.10) t i>[u]-£ vj[Fj(u)]=0 

>=i 

Here [•] denotes the difference between values in the brackets on the two sides of T. In the context 
of the above spectral Chebyshev approximation, the above equations read as follows: 


(2.11) t v(ujf u iW u %) ~ = 0 

;=i 

At each collocation point Ep on T (2.11) yields p equations for the 2p + 1 unknowns: ujy,uj y 
and tv. We need therefore p + 1 further independent conditions which will be provided by the 
compatibility equations, which, in the current situations, can be determined as follows. Let us 
define the characteristic matrix for ftj at the point P € T: 


( 2 . 12 ) = v 1 a Mn) 

y=i 

and denote by A* > and the eigenvalues and left eigenvectors, respectively, of A( i/ 1 ). Similarly, 

let us denote by A* 2) and t; (2) respectively the eigenvalues and left eigenvectors of the matrix 


(2.13) A(v 2 ) = '52(-vj)A j (u 2 n ) 

i = i 

(note that this time A*,) ^ -A* 2) and v* X) j v* 2) in general, since ti]y and xt 2 N are not coincident on T 
anymore). Assume that the (real) eigenvalues are ordered as follows 

A(l) < - ^(l) 1 > ^(2) < ^(2) > = 1, ...,p — 1 

We assume that T is the surface of propagation of a k-shock . so that the following entropy conditions 
are verified (e.g.[S] p.261): there exists an integer k e {1, ...,p} such that: 

Aq ) 1 <w< A*]) , A ( ‘ 2) <w< A ^ 1 


(2.14) 
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Fig. 2. 1 A k-shock with p=3 and k=l 

Figure 2.1 illustrates an example of a k-shock for a one-dimensional problem. For a fixed time t 
we have drawn in bold the straight line with slope = w (the speed of the shock front). On its left 
(resp. wright) we have labelled with j the straight line whose slope is equal to the characteristic 
speed Aj !) (resp. for / = 1, ...,p. 

The compatibility equations for Oj will be therefore given by the p - k + 1 equations 


*»> ' (l^ + A <» S^r) " ”«> ' _ ^ A ‘^' n) °dT r/ j 


(2.15) 

while those for Q 2 are given by the k equations: 


r = k,...,p 


(2.16) 


■( du "+\ t . f/_ Y'^r u 2)^L r .) 

(2) \~dT A(2) ^J/ ~ (2) I 7 dr J ) 


a =p- k+ l,...,p 


As usual, t is the unit vector tangential to T at the point P under consideration. 
Notice that (2.14) can be rewritten in the form (see [L], p.25) 


(2.17) Afo 1 <w<\$ , Af 2) <w< Af 1} 

which shows that there exists only one index k such that the shock speed w is intermediate to the 
characteristic speeds A* on both sides of the shock. 

In the case of a flow regime supersonic in fii (upstream subdomain) and subsonic in Q 2 (down- 
stream subdomain) (see, e.g., Fig.2.1) the flow field within is uninfluenced by that within £2 2 , 
and it should therefore be computed firstly. 
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Remark 2.2 (Discontinuity due to unsmooth initial data) 

For linear hyperbolic systems ( Aj independent of u in (1.1)) with discontinuous initial data, there is 
propagation of discontinuity across fronts whose position and motion can be a priori determined in 
terms of the data. In these cases, w is given and A* > = — A**), v*) = for all k = 1, ...,p, so that 
the p resulting compatibility equations together with the p Rankine-Hugoniot conditions (2.11) 
allow the calculation of the 2p unknowns ujy, ujq at each collocation point on the discontinuity 
fronts.o 

3. An Iteration-By-Subdomain Algorithm for the Solution of the Domain Decomposition 
Problem 

Let us go back to the case of solutions which are continuous across the subdomain interface T. 
The description of the domain decomposition method given in the previous section shows that the 
spectral collocation problems in fti and Q 2 are coupled throughout the continuity conditions (2.4) 
at the interface T. To remove this coupling we propose the following iterative method. Assume the 
solution is available at the n-th step. Let v, r, A* and v k be defined as at the point (c) of section 2 
for each collocation point of £r (clearly, A* and v k depend on the value of the solution at the points 
ofT). 

Let A (u) and T(u) denote respectively the matrix of the eigenvalues of A(u) and that of the corre- 
sponding (left) eigenvectors, so that: 


(3.1) T(v) A(u) = A (u) T(v) 

Then define at the point under consideration 

(3.2) X 1 = ( w lv)fr X 2 = ( w tf )|r 

Finally, we denote by T q (u) the q x p matrix given by the first q rows of T(y), while T p _,( v) will 
denote the (p — q) x q matrix given by the remaining p — q rows of T(y). The solutions at the new 
step, say (u]y ) n+1 in Qj and (ujy ) n+1 in ft 2 , can be obtained by solving two independent problems 
in fti and in ft 2 , respectively. 

Precisely, (uj v ) n+1 satisfies the interior equations (2.2), the interface equations (2.5) together with 


(3.3) Tp- q (v) (uir) n+1 = Tp- q (u) X 2 on T 

and, finally, the boundary equations given at the point (d) of sect.2. 

Similarly, (u^) n+1 satisfies the interior equations (2.3), the interface equations (2.6) together with 


(3.4) 


T q (v) (u^) n+1 = T q (u) x 1 onr 
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and, finally, the boundary equations prescribed at the point (e) of sect.2. Note that at the limit of 
the convergence process, the (independent) conditions (3.3) and (3.4) will ensure the fulfillement 
of the continuity requirement (2.4). 

Remark 3.1 The extension of the above method to decompositions with several subdomains is 
straightforward. At each step, one is left to solve with as many independent subproblems as the 
number of subdomains. All these subproblems can be solved simultaneously within a parallel 
computer environment. 

4. A Particular Case: One Dimensional Problems 

We consider the special case in which Q is the one-dimensional interval (-1,1), and fii = (— 1, a), 
Q 2 = (a, 1) for some 0 < a < 1. In this case T = {a}. 

The differential system reads as follows 


(4.1) ^ + A(u) — = / in O x (0, T) 

Denoting with {A*, k = 1, ...,p} and {v*, k = 1, ...,p} the eigenvalues and the left eigenvectors of 
A(u), respectively, the compatibility equations (1.7) in the current case become 

< 4 - 2 > t=i * 

Setting A = diag {A 1 , ..., A p } and denoting by T the matrix whose k-th row contains v k , we can 
write (4.2) as follows 


(4.3) T ^ +AT ^ = r/ 

Obviously T and A depend on u if so does A. 

Assume that at the interface point x = a the first q eigenvalues of A are positive. (Of course, 
A, A, T and q depend on the time level <). Let us denote by A, and T q the q x p matrices obtained 
suppressing the last p — q rows of A and T, respectively. Then the following q equations 


du du 

(4.4) Tq — +A q T — =T q f at x = a 

provide the q compatibility equations for fii. Similarly, denoting with A p _ g and T p _, the lower 
part of the matrix A and T respectively, the equations 
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(4.5) T p - q — + Ap-j T — = T p . q f at x = a 

provide the p - q compatibility equations for fi 2 - The compatibility equations for Q 2 at the right 
hand boundary point x = 1 are obtained in a similar way, and take the form (4.4), while at the left 
hand boundary x = —1 the compatibility equations for Qi take the form (4.5). 

The hyperbolic system (4.1) needs to be completed by the initial condition 


(4.6) u(x, 0) = <p(x) iGfl 

and by a set of boundary conditions ( q equations at the point x = -1, and p - q at the point x = 1) 
that we assume having the form 


B\u = g x at x = — 1 

(4.7) n 2 , , 

B 2 u = g atx = 1 

where B\ is a q x p matrix, B 2 is a (p - q) x p matrix, and g x and g 2 are two given vector functions. 
Remark 4.1 Since A T = TA from (3.1), instead of (4.3) we can write 


(4.3)' 


+ A 


du 

dx 


= Tf 


o 


This is exactly the form that was used in [CQ] to derive the numerical boundary conditions at the 
physical boundary of the domain. 

We also notice that if A is a constant matrix, then setting z = Tu (characteristic variables), from 
(4.3) we obtain the following characteristic form of the compatibility equations 


(4.3)" 


dz dz 
dt + ^ dx 


-T f 


o 


We will now apply the multidomain spectral collocation method described in section two to the 
one dimensional problem (4.1). Let us define 


= -l + 


a + 1 


(tj + 1) 


xj = a + 


~~ 2 ~~ + 1) > / = 0, ..., N where tj = cos 


2 


2 
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The points tj are the Chebyshev-Lobatto nodes in the reference interval [-1,1], while x] and xj are 
their images within the subdomains Qi and O 2 respectively (note that - 1 < x] ■ < a < xj < 1 for 
1 < j < N - 1). 

At each time t > 0 we look for tt}y(f) e (Pn(^i)) p > u^(t) € (Pn(^ 2 )) p satisfying: 

(a) internal equations 

(4.8) ^ + A(u),)^=/ atx) \<j<N-\ 


(4.9) d ft +A(u " )d f N =/ i < > < JV - i ; 

(b) interface equations 


(4.10) 

T ' at +A,T dx = T - 1 otl = “ 

(4.11) 

rp & U N . * rp & U N _rp f . _ 

1 p -9 Q t + A p-« 1 q x 2 p-i J at x a 

(4.12) 

«]yr = tijf at x = a ; 

(c) boundary equations 


(4.13) 

T ' d m +A ' Td l N=T ' f atx = l 

(4.14) 

rp & U N _ rp f at X ~ 1 

J P-« Q t + A P-« J - J P-iJ aZX 1 

(4.15) 

Bixtff =g 2 at x = 1 
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(4.16) B\u x N = g' at x = —l 

Note that (4. 10) and (4.14) are the compatibility equations for Qj at the points x = a and x = —l, 
respectively, while (4.11) and (4.13) are those of Sl 2 at the points x = a and x = 1, respectively. 
Finally, (4.15) and (4.16) are the boundary conditions (see (4.7)). 

Remark 4.2 (Fully discrete approximation) 

A fully discrete approximation to the problem (4. 1 ), (4.6), (4.7) can be obtained using a time march- 
ing scheme in (4.8)-(4.16). Whatever scheme (either implicit or explicit) one adopts to advance 
from a known time level t k to a new one t fe+1 , the matching interface condition (4.12), as well as 
the boundary conditions (4.15) and (4.16) should be enforced at the new time level. 

If an explicit time stepping is used, at the time level t k+1 the unknown vectors «(*]•)} and 
{u^(xy)},7 = 1 , ..., N - 1, can be computed using solely the internal equations (4.8) and (4.9), re- 
spectively. Once these internal values are available, the interface equations (4.10)-(4.12), together 
with the boundary equations (4.13)-(4.16), can be solved to provide the remaining values {u] v (x] )} 
and (u^(x?)} for »' = 0 and N. Actually, we note that the presence of derivatives in space among 
boundary and interface equations relates boundary and interface values to each other. We also em- 
phasize that the differential equations (4.10), (4.11), (4.13) and (4.14) ought to be advanced by the 
same explicit time stepping that was used for the equations at the internal points. 

When an implicit time stepping is used, the internal unknowns are not decoupled from the remaining 
ones anymore. We will see an example in the next section.^ 

We write now the iteration-by-subdomain method described in section 3 for the solution to the 
problem (4.8)-(4.16). Assume (t»Jy) n , (u^) n are available at the n-th step and define 


(4.17) X 1= (t*]*r) n , X 2 = («*)" at x = a 
Then in Qj we look for (ujy) n+, (f) G (Pn(Qi)) p such that 

(4.18) ~ (xj,)"*' + A ^ <«] = / atx) , 

(4.19) T,y t («),)"*' +A,r^(«i,r 1 =r,/ at x = a 


(4.20) 


T p - q (t iJr) n+1 = T p . q x 2 at x = a 
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(4.21) v, ^ (u'„r' + Ap-, t A (u‘ K r' = t,-, i (it % ~ i 


(4.22) B, («},)"*'= g' atx = - 1 

where the matrices A, A and T depend on (ujyr) n+1 . 

In £^2 we solve for (u^) n+1 (t) G (Pn(£22)) p satisfying 

(4.23) ^ ( 4 )”* ' («Jr)" +1 - / «< , 1 < j < N - 1 


(4.24) ^ (uj,)” +1 +A,-, (uj,)”* 1 =T^, / atx = a 


(4.25) T g (u^) n+1 =T q x X at x = a 


(4.26) T„ ^ (u^) n+1 + A t T^ (4)" + 1 = T, / a* x = 1 


(4.27) B 2 (u^) n+1 = p 2 ot x = 1 

In (4.23)-(4.26) the matrices >1, A and T depend now on («^) n+1 . 

Note that the problem in Q] is independent of that in fl*- The extension of the above iteration-by- 
subdomain method to the case of a subdivision by M subdomains ( M > 2) is straightforward. At 
each iteration we obtain M independent subproblems that can be solved simultaneously. 

Here above, the iteration-by-subdomain method has been applied to the semidiscrete (continu- 
ous in time) problem (4.8)-(4. 16). To be effective, the iteration method should be applied at each 
time-step after that a time marching scheme has been used to get a full space-time discretization 
of the problem (see the previous remark). In this way, at the new time-level t* +1 , we can ap- 
ply the iterative method until when we achieve convergence to the spectral multidomain solution 
u l N (t k+ '),u 2 N (t k + l ). 
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5. Convergence Analysis for the Iteration by Subdomain Method 

In this section we will prove that the iteration by subdomain method introduced so far is convergent. 
Our analysis is confined to one-dimensional hyperbolic systems with constant matrix. The initial- 
boundary value problem under investigation takes the form (4.1), (4.6), (4.7) where A is a 2 x 2 
matrix which is not restrictive to assume of the following form 


(5.1) 



with |a| < 1 


The eigenvalues of A, say A and —fi, have opposite sign, as A = a + 1 >0 and —fi = a — 1 < 0. 
We notice that in this case 


(5.2) 


A = T AT , where A = diag {A, -fi} and T = -^= ^ j ^ 


Among the various sets of boundary conditions that render this problem well posed we consider 
the following ones 


(5.3) «i(-l,t) = 0 , «i(M) = 0 Vt€(0,T) 

We will require the initial value <p( x) (see (1.2)) to be a continuous vector function, with a first 
component vanishing at x = — 1 and x = 1 . Under this assumption, initial and boundary conditions 
are compatible, hence the solution to our problem is continuous for all time. Assume that we are 
advancing in time the spectral multidomain problem (4.8)-(4.16) by an implicit method. For the 
sake of simplicity, let us consider the backward Euler method (no essential modification occur 
where considering another implicit method). When we advance from t k to t k+l , at the new level 
t k+] the new functions w 1 := («Jy)(<* +1 ) and u 2 := (n^ r )(f fc+1 ) satisfy the following multidomain 
problem: 


(5.4) 


pu 1 +Au l x = f^ rev 


at Xj , 1 < ; < N - 1 


(5.5) 




at x 2 , 1 < j < N — 1 


r, /? u 1 + A r 9 = T q fan 


(5.6) 


at x = a 
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(5.7) T p _, p u 2 - n T p _, u\ = T p _, at x = a 

(5.8) T q u 2 = T q u 1 at x = a 

(5.9) r p _, u 1 = r p _, u 2 x = a 

(5.10) T, /? u 2 + A T, u 2 = T q /* w atx = l 

(5.11) T p - q ft v} — nT p - q u\ = T p - q fprg,, at x = — 1 


(5.12) tt{ =0 at x = — 1 


(5.13) «? = 0 at x = 1 

In the previous equations p is the inverse of the time-step, where denotes 

the value of ujy at the former level t k , and f 2 rev is defined similarly. Note that in the current case 
T q and T p - q denote the first and second row of the matrix T respectively. The two equations (5.8) 
and (5.9) correspond exactly to the continuity condition (4.12). 

Let us apply now to the problem (5.4)-(5.13) the iteration-by subdomain method described in the 
previous section. At the (n+l)-step, the solutions (u 1 ) n+1 and (u 2 ) n+1 satisfy the equations (5.4)- 
(5.7) and (5.10)-(5.13); in accordance with (4.20) and (4.25) the matching interface conditions (5.8) 
and (5.9) are dealt with as follows 


(5.14) 

t 


\ 

i 

i 


T,(u 2 ) n+1 =r,(u‘) n at x- at 


(5.15) 


r p _,(ti , ) B+1 = Tp. q (u 2 ) n at x - a 
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The main result we are going to prove in this section is that the above iterative procedure converges, 
i.e., (u 1 ) B — ► tt 1 and (u 2 ) n — ► u 2 as n — ► oo. Furthermore, the convergence rate is independent of 
the polynomial degree N. 

The proof is rather involved; it is convenient carrying out it using the characteristic form of the 
equations. For that, let us define the error functions (in characteristic form) at the (n + 1) — th 
iteration as follows 


(5. 1 6) (e‘) n+1 := T[(u‘) n+1 - «’] for i= 1,2 

By comparing the equations (5.4)-(5.13) with the equations satisfied by the iterates (u‘) n+1 , t = 1, 2, 
it is readily seen that the following set of error equations hold: 

within f2] (as usual, subindex denotes component): 

(5.17.1) /?(e’) n+1 + A (e 1 )* +1 =0 at x) , 1 < j < N - 1 

(5.17.2) (e j ) n+1 + (4) n+1 =0 at x = -1 

(5. 1 7.3) /?(e!) n+1 - M(« 2 )* +1 =0 at x = - 1 

(5.17.4) /?(ei) n+1 + A(e|)" +I =0 at x = a 

(5.17.5) ( C i)" + i = ( C 2) n atx = cx 
within Q 2 : 

(5. 18.1) /9(e 2 ) n+1 + A(e 2 )” +1 =0 at x) , 1 < j < N - 1 


( 5 . 18 . 2 ) 


(e 2 ) n+1 + (e 2 ) n+ ' =0 at x = 1 
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(5. 1 8.3) /?(e?) n+1 + A(e?£ +1 =0 at x = 1 

(5.18.4) P(e 2 2 ) n+l - /*(ci)” +1 =0 at x = a 

(5. 1 8.5) (c?)” +1 = (e{ ) n at x = a 
We start proving two important lemmas. 

Lemma 5.1 Let tj > 0 and e be two given real numbers, and consider the following "backward" 
collocation problem in the reference interval [-1,1]: find u € Pn such that 

(5.19.1) u — i;u x =0 at tj,j=l,...,N 


(5.19.2) u = e at to = 1 

where, as usual, tj = cos%,j = 0, ..., N are the Chebyshev-Lobatto points in [-1,1]. There exists 
a constant a ffiv) whose absolute value is strictly less than 1 such that 


(5.20) u(— 1 ) = a N (ij)e 

Proof Let Tk(x ) denote the Chebyshev polynomial of degree k. (We recall that T*(x) = cos k9, 
where $ = arcos x). Since for j = 1, ..., N — 1, tj are the roots of T' N (x), from (5.21 . 1) we deduce 
the following identity: 


(5.21) u - ij u s = r(^j\r_i + <f> N ) VxG[-l,l] 

where i(x) = T' N (x), 4>n(x) = xT N (x), and t is a constant that can be determined using the 
inflow boundary condition (5.19.2). Following [C] we can state that 

(5.22) degree (4> n ) = n, parity of <f> n = parity of n , (fi n )k > 0 

for n = N - 1 ,N and k = 0, n. Here (4> n )k are the Chebyshev coefficients of <f> n , i.e.. 
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(5.23) (0, 

Owing to (5.21) we can set 


2 f x 1 

’„)* = / <f> n (x)T k (x) — dx with c 0 = 2 and cjt = 1 if fc > 1 

tfc* J_i Vl - x 2 


(5.24) « = r(<p v + V»,) 

where <p v and tp,, are the unique solutions of the following ordinary differential equations: 

(5.25) ^ € Pn : <Pti ~V ,% = ^n-i 


(5.26) V'ij ^ Pn : V>«? ~ V V’ij.x = <£n 

If p is any polynomial of Pn we recall that (e.g., [CHQZ] Ch.2): 


N 


JV-1 


(5.27) p(x) = V pjt r fc (x) , p x (x) = J2 Pk l) T k(x) tip = ~ Y) m Pm 

fc=0 *=o c * ->*♦• 


where the symbol m k means that |m — k\ is odd. Taking into account (5.27), from (5.25) and 
(5.26) we obtain the following recurrence relations for the Chebyshev coefficients of ipr, and W 


(5.28) 


(<P V )m = (^AT-l)m + V X) W* m = N,N — 1, 0 


k>m+l 


(5.29) (ij> n ) m = (^Jv)m + *1 “ S m = N, N — 1, 0 

Owing to (5.22) we conclude that 


A^m+1 

tq/lm 


(5.30) (<pr,) m > 0 , (tp n ) m >0 for m = 0, ..., JV , for all r}> 0 

Moreover, noting that T*(l) = 1 for all k from (5.24) we obtain 
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N 

(5.31) r = sf + ($ v ) m ) whence sign r = sign e 

m=G 

It follows that the solution to the collocation problem (5.19) has the following Chebyshev coeffi- 
cients 

N 

(5.32) u m = ((£,) m + (^,) m ) e/ Y, ((£*»)"* + (^)m) 

m=0 

and therefore 


(5.33) sign = constant = sign s , m = 0, ...,N 

Noticing that ti(-l) = £iU 6 m(-l) m as T m (-l) = (-l) m , we deduce (5.20) from (5.32) by 
setting 

N AT 

(5.34) ff ff (t ;) := £(- 1) TO (( <p,) m + (fej/ Y + AU 

w=0 jfe=0 

Finally, the property 


(5.35) KM<1 Vt; > 0 

follows from (5.30). o 

Lemma 5.2 Lef rj > 0 and e be two given real numbers, and consider the following "foreword" 
collocation problem in the reference interval [-1,1]: find vcPn such that 


(5.36.1) 


v + rjv x = 0 attj , j = 0, ...,JV — 1 


(5.36.2) 


v = e attff = 


where the points tj are defined as in Lemma 5.1. Then 
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(5.37) v(l) = a N (r))e 

where on w the same as in (5.34). 

Proof We notice first of all that 


(5.38) v + ijv s = p(to-i - to) Vx G [- 1, 1] , p = 

The constant p has been determined by the boundary condition (5.36.2) noting that 7^(1) = k 2 and 
T' k (—\) = (— l) k+i k 2 . Consider now the following "backward" problem associated with (5.36): 
look for u; G Pn s.t. 


v)-rjw x = 0 at tj , j = 1, JV 
w = e at to = 1 


We want to show that 


(5.39) v(— x) = w(x) Vx G [—1, 1] 

As a matter of fact w verifies: 


(5.40) w - i) w x = p*(<t>N-\ + <!>n) Vx g [—1, 1] , p* = - — ~ 

Writing (5.40) at the point — x we obtain 

tu(-x) - 17 w x (-x) = p\ 1 - x) r N (~x) Vx G [-1, 1] 
whence, in view of (5.39), 


(5.41) v(x) + r, v x (x) = [(-1)^- Vl (1 - x)r N (x) 

We have used the property that 1*(— x) = (-l)* +, l*(x) Vx G [—1, 1] and the obvious relation: 
w x (-x) = - v,(x). Since (- 1)^ - V* = P and (1 ~ x)Tff(x) = 4 >n-\ - toy we conclude from (5.41) 
that v satisfies (5.38), whence v is the only solution of (5.36). 

We can now apply the previous lemma to get tv(— 1) = a ff(ij)e and therefore (5.37) follows owing 
to (5.39). 
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Furthermore, since the Chebyshev coefficients of w are given by (5.32), using again (5.39) we 
deduce that the Chebyshev coefficients of v are 

AT 

(5.42) V m = (- l) m [(£,) m + (iMml e/ X) ((^)m + A)rn) 

m=0 

In view of (5.20) and (5.37), the factor defined in (5.34) will be called the outflow! inflow ratio.o 
We are now in the position to state the following result 
Theorem 5.1 If we define 

(5.43) A' = 2 \/P( 1 + a) , p! = 2/i//3( 1 + a) 
f/ie solution to the problem (5.17) satisfies : 

(5.44) (e{ ) n+1 (a) = -ajv(A') <r*(/i') ( e\) n (a ) 

Proof Let t(x) = ^j(x + 1) — 1 be the affine transformation from Qj = [— 1, a] to the reference 
interval [-1,1]. Let us define: 

(5.45) « G P N : u(t(x)) = (el) n+1 (x) Vx £ Q, 

It is readily seen from (5.17.1), (5. 17.3) and (5.17.5) that u is the solution to a backward collocation 
problem like (5.19), provided we set tj = p! and e = (e^Tia). By (5.20) we therefore deduce 

(5.46) (e!) n+, (-l) = ti(-l) = a N (yl) (e\) n (a) 

Let us define now: 


(5.47) v G P N : v(t(x)) = (ej ) n+1 (x) Vx £ 

Owing to (5.17.1), (5.17.4) and (5.17.2), it follows that v is the solution in [-1,1] of the foreward 
problem (5.36), provided we set rj = A' and e = — (e 2 ) n+1 (— 1). We can therefore apply the result 
(5.37), and owing to (5.47) and (5.46) we obtain 

(«1 ) n+1 (or) = v(l) = <rjv(A')e = -<ryv(A')<r^(/i')(c 2 ) n (a) 
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The following result is the counterpart of theorem 5.1 for the subdomain Cl 2 . 

Theorem 5.2 Let us set 

(5.48) A" = 2A//?(1 -a) and p" = 2p/P(l - a) 

The solution to the problem (5.18) satisfies 

(5.49) (^r'(a) = -<r A ,(A") »*(/.") (e|)”(a) 

Proof Let t(x) = - a) - 1 be the affine mapping from Q 2 = [a, 1] into [-1,1]. and define: 

(5.50) v G P N : v(t(x)) = (e?) n+1 (x) Vx G U 2 

Owing to (5.18.1), (5.18.3) and (5.18.5) we deduce that in [-1,1] v satisfies a problem like (5.36) 
with tj replaced by A" and e by (e{) n (a). In view of (5.37) we have therefore 

(5.51) (cf ) n+1 (1) = ffjv(A") (e]) n (at) 

On the other hand, from (5.18.1), (5.18.2) and (5.18.4) it follows that the function 

(5.52) u G P N : u(t(x)) = (el) n+1 (x) Vx G U 2 

satisfies upon [-1,1] a problem like (5.19) provided ij is replaced by p" and e by — (e,) n+1 (l). Thus 
from (5.20) we obtain 

(el) n+1 (a) = -^( M ,, )(c?) n+1 ( 1) 

Now (5.49) follows from (5.51).o 

Let us define the following sequence of interface errors: 

(5.53) E n = [((e{ ) n (a)] 2 + [(ei) n (a)] 2 + [(e?)" (a)] 2 + [(t\) n (a)] 2 for n> 1 
From the previous theorems we deduce the following convergence result. 

Theorem 5.3 The interface error defined by (553) reduces at each iteration according to the law 
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(5.54) £ n+1 < o' N ( A, /*; a)£ n , n > 2 
w/zere r/ie reduction factor is defined as follows 

(5.55) »J,(. A, /i; a ) iaVjUA ^(A-yJrOO} < 1 

Proof Owing to (5. 1 7.5), (5. 1 8.5) and to the theorems 5.3 and 5.4 we obtain that the interface errors 
propagate according to the following relation 


(5.56) ET +X ={[(ei) n ] 2 + [(e^) rt ] 2 }(r^(A>^(^ + {[(e!) n ] 2 + [(ci) n ] 2 }ffAr(A>^(/i") forn>2 


The inequality (5.54) follows easily. Note that <r^(A,/i;a) < 1 for all positive A and p as a 
consequence of (5.35). ❖ 


Remark 5.1 (Behaviour of the error reduction factor). 

The behaviour of o* N ( A, fx; a) (and, by consequence, that of the interface error sequence (5.53)) is 
driven by the behaviour of the outflow/inflow ratio defined in (5.34). 

From tables 5.1 and 5.2 below we see that for all values of N and rj, a n(v) > 0* Moreover for any 
fixed value of t],0N(v) is uniformly bounded from above by a constant strictly less than one as N 
increases (even better, the value of o is substantially independent of N). This property is very 
important for it ensures that the error reduction factor does not approach one as N tends to infinity, 
yielding therefore a convergence rate for our iterative procedure which is (practically) independent 
of N. In table 5.2 we report the limit of Offiv) as N tends to infinity, for several values of rj. It is 
apparent that 


(5.57) 


lim on( t?) = e V rj > 0 

n/oo 


Thus, the limit is precisely the outflow/inflow ratio of the differential case, i.e., the value 
v(- l)/u(l), where v is the solution to the ordinary differential equation v — rjv s = 0. 

On the other hand, for fixed N, <rjv(rj) behaves like a monotonically increasing function of 17 (if 
N and/or 17 are not too small). In the current application, 17 takes the values of A', A", (f and fx", 
which are all proportional to the time-step At = /9 _1 (see (5.43) and (5.48)). Therefore, for large 
N we have approximately 


<r^(A, 71 ; a) ~ exp{-^2a*(A + fx)/Xfx} 


a* = min ( 1 + a, 1 — a) 


(5.58) 


26 


We recall that A and -/i are the eigenvalues of the transport matrix A (see (5.1)), while x = a is the 
abscissa of the interface between the two subdomains £2i and Qi- Notice that a* is the minimum 
measure of the subdomains.o 


X 

4 

8 

16 

32 

0.01 

0.05 

0.1 

0. 5 

1 . 

5. 

10. 

100. 

0.10195409 

7.6327433E-2 

5.8272078E-2 

2.8037383E-2 

0.13687601 

0.67032260 

0.81873085 

0.98019867 

4.532901 8E-2 

1.8829854E-2 

6.00765 84E- 3 

1.8321630E-2 

0.13533534 

0.67032005 

0.81873075 

0.98019867 

2.5022730E-4 

1.975305 IE-4 

1.7829169E-6 

1.8315638E-2 

0.13533528 

0.67032005 

0.81873075 

0.98019867 

9.4749803E-11 

1.8740999E-11 

2.061 1536E-9 

1.8315638E-2 

0.13533528 

0.67032005 

0.81873075 

0.98019867 


Table 5.1 The value of the outflow/inflow ratio (n)for several values of rjand N 


T1 

<L Cri) 

0.01 

0. 

0.05 

4.24835425 E-18 

0.1 

2.0611536 E-9 

0.5 

1.8315638 E-2 

1. 

0.13533528 

5. 

0.67032005 

10. 

0.81873075 

100. 

0.98019867 


Table 5.2 The value of cl 


(T)) := lim (-p) for several values of i) 


From the previous theorem we have that 


(5.59) /im n _ 00 .E' n = 0 

In order to conclude that our iteration-by-subdomain procedure converges it remains to prove that 
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the error functions (e*) n (x), * = 1,2, defined in (5.16) attain their maximum values at the interface 
point x = a. This is stated by the following theorem. 

Theorem 5.4 With the notation (5.16) and (5.23) we have 


(5.60) 


mox_,<*< a {[(e{ ) n (x)] 2 + [(e}) n (x)] 2 }+ 
max a < x <\ {[(e?)"(x)] 2 +[(e 2 ) n (x)] 2 } < E" V n > 1 


Proof From the proof of theorem 5.1 we easily deduce that for all x e [- 1, a]: 


N 


(5.61) (ei) n+, (x) = ^ti m T m (t(x)) = 


(e 2 )"(<*) 


N 




m=0 




and 


(5.62) 

(e{) n+, (x) = X} T m «(x)) = -<r *(#*') 

171=0 


(e|)"(«) 

EfooK0v)» + (l5v)t] 


N 

Y^UMrn + (i>\')m]T m (t(x)) 
m= 0 


Owing to (5.30) and the fact that |T m (£)| < 1 G [-1, 1], we deduce from the previous equalities 
and from (5.17.5) that for all x € [— 1, a] 


(5.63) |(ci) n+, (x)| < |(4) n+, (a)| , |(e!) n+, (s)| < M/*')||(e2)" +, («)| 

In a similar way, using theorem 5.2 we prove that for all i€[«,l] 


(5.64) |(e?) n+, (*)| < |(c?) n+, (a)| , \(e 2 2 ) n+ \x)\ < MA")||(e?) n+1 (a)| 

The inequality (5.60) follows easily from (5.61)-(5.64).o 
We can now state our final convergence result. 

Corollary 5.1 The iteration-by-subdomain procedure applied to the multidomain Chebyshev col- 
location problem (5.4)-(5.13) converges as n -* oo. 

Proof The result follows from (5.60) and (5.59). Actually, since T is non singular, the convergence 
of the sequence (e*) n implies that of (ti’) n — u * owing to (5.16).o 
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Remark 5.2 The same kind of result holds if the spectral Legendre collocation method is used 
instead of the Chebyshev one. 

Remark 5.3 (Convergence for decomposition with several subdomains) 

For a decomposition of [-1,1] using more than two subdomains we can essentially carry out the 
same type of convergence proof. Aside the technical difficulties, it is not hard to see that in this 
case the error reduction factor at each interface behaves like (5.58), where now a* is the measure 
of the smallest subdomain of the decomposition. 

If for instance, the computational domain [-1,1] is partitioned into M subdomains of equal length, 
the error reduction factor behaves like 


f 1 4 X+n 1 

(565 > 

hence the convergence rate slows down as far as M increases.^ 

6. Capacitance matrix interpretation 

We construct now the capacitance (or influence) matrix associated with the multidomain problem 
(5.4)-(5. 1 3). First of all we rewrite the problem in terms of the characteristic variables 

(6.1) a? 1 = Tti 1 , z 2 = Ttt 2 


If we define 

= Tf^ and h 2 = Tf from (5.4)-(5. 13) we obtain: 

(6.2.1) 

p z\ + A zj x = h\ at Xj , j = 0, ..., N — 1 

(6.2.2) 

pz\- fM z\ x = h\ at xj , j = 

(6.2.3) 

z \ + z \ - 0 at x\f = — 1 

(6.2.4) 

= z\ at x = a 


(6.2.5) 


at x = a 



( 6 . 2 . 6 ) 


z\ + z \ = 0 at Xq = 1 


(6.2.7) P z\ + A = h \ at x) , j = 0, ..., J\7 - 1 


(6.2.8) pz\- nz\ x = h\ atxj , j = 

Note that (6.2.4) and (6.2.5) are equivalent to (5.9) and (5.8), respectively. The capacitance (or 
Schur complement) system is the one that allows the determination of the values of the two char- 
acteristic variables z\ and z\ (and henceforth, those of z\ and z\ owing to the continuity relations 
(6.2.4) and (6.2.5)) at the interface x = a. 

Consider the polynomial solutions to the two following multidomain collocation problems: 


I p U\ + A U\ x =0 at x) , j = 0, ..., N - 1 
pU\-tiU{ x =Q at Xj ,7 = 1, ...,N 
U\+U\= 0 atxl = - 1 
U\ = 1 at x = a 


{ pUf + X U* s = 0 atx], j = 

PU%-n Ul >x =0 at x) , j = 1 , ...,N 
17 , 2 + Ui= 0 at xl = 1 

Uf = 1 at x = a 

Furthermore, let us consider the two polynomial solutions to the multidomain collocation problems: 


{ PV;+XV,\ x = h\ at xj , j = 0, ...» iV — 1 
PVi -t*V 2 \ x = h\ at Xj , j = 1,..., JN7 
V, 1 +VJ =0 at Xff = —\ 

Vj = Oat i = a 


’ /? Vj 2 + A V, 2 X = h\ at xj,j = 0, ... f N — 1 
, PVi-nVi tX = hl at xj , j = 
Vi 2 +y 2 2 =0 at x = a 
V] 2 = 0 at x — at 


( 6 . 6 ) 
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Let us denote by £ the (unknown) value of z\ (= z\) at x = a, and by rj the (unknown) value of 
z\ (= z%) at x = a. The following relationships hold between the solution to the problem (6.2) and 
those of the auxiliary problems (6.3)-(6.6): 

(6.7) z 1 = V x + v U x in Qi , 


(6.8) z 2 = V 2 + £U 2 inQ 2 

The values of £ and rj are therefore determined by the two following equations that arise from 
(6.2.4) and (6.2.5): 


(6.9) 


V, 1 (a) + r? U{ (a) = V 2 (a) + £ Ufa) , 


(6.10) Vfa) + v Ufa) = Vfa) + e Ufa) 

Using the last equations of (6.3)-(6.6) we obtain the 2 x 2 system: 


( 6 . 11 ) 



The matrix in (6.11) is the capacitance matrix (it is also called influence or Schur complement ma- 
trix) ; hereafter it will be denoted by S. 

In order to give a precise form to its entries, we note that by comparison of (5.17) with (6.3) U x 
and (e 1 ) n+1 satisfy the same set of equations in Qi. The only difference regards the value attained 
at x - a by the second component of the solution (precisely, U\(ci) = 1 and (e 2 ) n+1 (a) = (e 2 ) n (a)). 
Owing to (5.44) we can therefore conclude that 


(6.12) -U}(a) = o N (\')<r N (ti') 

with A' and y! given in (5.43). 

Exploiting a similar analogy between U 2 and the solution of problem (5.18) we easily obtain from 
(5.49) that 


(6.13) 


-UiM^anmapriy”) 



31 


where A" and p" are defined in (5.48). Then 


(6.14) 


S = 


( ' 


ffw(A>jv(/x') 

1 


) 


Theorem 6.1 The capacitance matrix is positive definite. Moreover, if a = 0 (i.e., and Q 2 have 

the same measure), then S is symmetric and its eigenvalues are 


(6.15) \ h2 = \±o N (\ , )<r N (p') 

Proof S is positive definite since |<r^(rj)| < 1 for any positive rj. The eigenvalues of S are 


(6.16) A , , 2 = 1 ± ^Jv(A') W) ' NiK) o N (p”) 

If meas fii = meas CI 2 (in our example, this means that a = 0 ), we have A' = A"(= 2A/0) and 
p 1 = p"{= 2 /x//?), whence 


frjv(A') = <?n( A") , on(h') = on(p") 

Thus S is symmetric, and its eigenvalues are given by (6.15) owing to (6.16).o 

We conclude this section by establishing an important relationship between the matrix S and the 
iteration-by-subdomain method we described in the previous section. 

Let us consider the spectral multidomain problem (5.4)-(5.13), and denote by 5 = (£,* 7 )* the un- 
known vector of the characteristic solution at the interface point x = a, i.e.. 


(6.17) £ = z\ (a) = z?(a) , ij = z\{a ) = z\(a) 

By a straightforward calculation it is not hard to see that one step of our iteration-by-subdomain 
procedure amounts to apply one step of the Richardson iterative method to the capacitance system 
( 6 . 1 1 ). Precisely, going from the step n to the step n+ 1 produces a change on the interface variables 
E which is ruled by the following relationship: 


(6.18) E n+1 =E n + (R-SE n ) , 

Here we have denoted by R the right hand side of (6.11). 

Using in (6. 1 8) an acceleration parameter w a priori different than 1 , i.e., considering the following 
sequence 
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(6.19) S” +1 = 5 n + u(R - S 5") , u > 0 

amounts to replace the interface matching conditions (5.14)-(5.15) by the new ones 


(6.20) T g (u 2 ) n+1 = uT q (u' ) n + (1 - w)T q (u 2 ) n a tx = a 


(6.21) Tp_ f (u , ) n+1 = u>T p _ 9 (u 2 ) n + (1 - w)T p _ g (u 1 ) n at x = a 

The other statements of the iterative algorithm hold unchanged. In (6.20)-(6.21) u plays the role of 
a relaxation factor. In view of (6.19), it is well known that whenever S is symmetric and positive 
definite, the best choice of u is (see [Y]) 


^opt — 2/(A m i n + A max) 

where A mtn and Xma X are the minimum and maximum eigenvalues of S. Thus if meas (Qi) = meas 
(Q 2 X in view of theorem 5.5 we conclude that u>opt = 1. 

7. Numerical Results 

We show the potential capability of our domain decomposition method to yield effective solutions 
of hyperbolic systems by presenting some numerical results for the simple problem considered in 
section 5. We allow however non homogeneous right hand side of (5.3). 

The time discretization method considered is the second order Crank-Nicolson scheme. 

We considered the case of two different sets of initial and boundary data, so that the corresponding 
exact solutions are: 


latest function: u 1 (x, t) = e * cos ax , u 2 (x, t ) = e 1 sin ax 


2 n< *test function: u’(x,t) = e t arctg (3(x + 0.5), u 2 (x,t) = e* arctg f3(x — 0.5) 
where a and /? are some given parameters. 

In all cases, we display the results obtained at the time t = 1 for several values of the time step A t, 
the polynomial degree N used inside each subdomain, and the number M of subdomains. In tables 
7.1 and 7.2 we report the maximum norm of the relative error between the approximate and exact 
solutions at the time t = 1 using At = 0.01. The tables refer to the first and second test functions 
respectively, for single domain approximations and for approximations with two subdomains of 
equal length. 
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one domain 



two subdomains 

A 

1 

4 

10 


N\ 

1 

4 

10 

4 

1.361E-4 

2.429E-2 

0.184 


B 

1.599E-5 

1.059E-2 

0.747 

8 

1.472E-7 

8.134E-3 

3.741E-1 



1.436E-7 

1.311E-5 

7.266E-2 

20 

1.472E-7 

1.546E-7 

3.718E-4 


R 

1.363E-7 

1.548E-7 

1.552E-3 


Table 7.1 First test function, one and two subdomain approximations 


X 

one domain 


X 

n\ 

two subdomains 

1 

4 

10 

1 

4 

10 

4 

4.767E-4 

1.288E-2 

2.812E-2 


B 

1.937E-4 

6.119E-3 

3.668E-2 

8 

1.094E-4 

5.079E-3 

2.345E-2 


H 

7.377E-7 

1.024E-3 

1.893E-2 

20 

2.572E-7 

2.860E-4 

1.005E-2 


H 

2.577E-7 

1.943E-4 

9.23 IE-3 


Table 7.2 Second test function, one and two subdomain approximations 

As we can see, the better accuracy we obtain using two subdomains rather than simply one is 
more relevant if the expected solution exhibits important oscillations (and/or variations) within the 
computational domain. Furthermore, since the matrices associated with the spectral collocation 
method are full, the use of several subdomains, accompanied with a sound tuning of the polynomial 
degree within each subdomain, allows the reduction of the overall complexity of the numerical 
problem. 

We report now the results obtained by the iterative procedure presented in section 3 for the resolu- 
tion of the multidomain problem. 

In the tables below we report the number of iterations that are needed in order to damp the initial 
error by a factor of 10~ 5 . As usual M denotes the number of subdomains and N the polynomial 
degree of the discrete solution inside each subdomain. 

Table 7.3 refers to the second test function when /? = 10: it shows that the number of iterations is 
independent of N, while it grows at most linearly with M when At is "large". According to the 
convergence theory of section 5, the rate of convergence is as larger as the time-step At is smaller. 
Finally, in the table 7.4 we present the results for the first test function when a - 1 , obtained for 
N = 16 and several values of M and A t. 
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\M 

N \ 

2 


8 

12 

4 

2(2) 

5(4) 

mm 

mm 

12 

2(2) 

4(4) 

cm 

B£Hfl 

20 

2(2) 

4(4) 

m 

m 


Table 7.3 Number of iterations for the second test function when (3 = 10 and At=0.1 
(within the brackets are the values for the case A t = 0.01) 



Table 7.4 
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